function [P,zeta] = dreumStationary_6_WTP(alpha,H,impossible,K,phi,s,tau,Vb)
                                           
Nutp=size(Vb{1},3);

% Approximate site choice probabilities for each applicant
P=cell(1,tau);
Pnot0 = cell(1,tau);

for idx2 = 1:tau

    PO = zeros(H+1,K,Nutp);

	for idx3 = 1:H+1
	
		for idx4 = 1:K
            
			if impossible(idx3,idx4) == 0
				                                     
				PO(idx3,idx4,:) = integral(@(z)Pint2(z,H,idx3,impossible(:,idx4),phi(:,idx4),...
                    s,Vb{idx2}(:,idx4,:),Nutp),0,Inf,'ArrayValued',true,'AbsTol',1e-12,'RelTol',0);
                               
            else
                
                x = (repmat(Vb{idx2}(1,idx4,:),[H,1,1]) ...
                    - Vb{idx2}(2:end,idx4,:))./repmat(phi(2:end,idx4),[1,1,Nutp]);            % Calculate preference shock
                
                CDF = [exp(-exp(-x(1:end-1,:,:))); 1 - exp(-s*x(end,:,:))];
                CDF(repmat(impossible(2:end,idx4),[1,1,Nutp])) = 1;
                PO(idx3,idx4,:) = prod(CDF);
                
 			end
		end
	end
 
% 	% Adjust approximated choice probs so the sum to 1 (Note: type 1
% 	% choice probs may sum to value > 1 due to approximation error.
% 	% This code normalizes the choice probs by the sum of all choice
% 	% probs such that they sum to exactly 1

    % Note that for people w/ max pp stock, Vb(0,K,:) = Vb(23,K,:) implying
    % the choice probability for j = 0 is zero
    z = PO;
    z(repmat(impossible,[1,1,Nutp])) = 0;
    z(1,:,:) = PO(1,:,:);
    z = z./(repmat(sum(z),[H+1,1,1]));
    y = repmat(z(1,:,:),[H+1,1,1]);
    z(repmat(impossible,[1,1,Nutp])) = y(repmat(impossible,[1,1,Nutp]));
    P{idx2} = max(z,1e-100);
%     P{idx2} = z;
    Pnot0{idx2} = z(2:H+1,:,:);
    Pnot0{idx2}(repmat(impossible(2:end,:),[1,1,Nutp])) = 0;

		
end
             
% Calculate stationary transition matrix A 
A = cell(1,tau);
for idx1 = 1:tau
    
    for idx2 = 1:Nutp
        
        a = zeros(K,K);
        
        for idx3 = 1:K
            
            a(idx3,min(idx3+1,K)) = P{idx1}(1,idx3,idx2) + (1 - phi(2:end-1,idx3)')...
                *Pnot0{idx1}(1:end-1,idx3,idx2) + P{idx1}(end,idx3,idx2)*(idx3==K);
            a(idx3,1) = phi(2:end-1,idx3)'*Pnot0{idx1}(1:end-1,idx3,idx2) ...
                + P{idx1}(end,idx3,idx2)*(idx3==1);            
            if idx3 > 1 && idx3 < K
                a(idx3,idx3) = P{idx1}(end,idx3,idx2); 
            end
                       
        end
        
        A{idx1}(:,:,idx2) = a;
        
    end
    
end


% Calculate stationary preference point distribution, zeta
zeta = cell(1,tau);
b = [alpha; zeros(K-1,1)];
for idx1 = 1:tau
    for idx2 = 1:Nutp
        zeta{idx1}(:,:,idx2) = inv(eye(K) - (1 - alpha)*A{idx1}(:,:,idx2)')*b;
    end
end


% kappa = cell(1,tau);
% for idx1 = 1:tau
%     for idx2 = 1:Nutp
%         kappa{idx1}(idx2,:) = zeta{idx1}(:,1,idx2)'*pit(idx1)...
%             ./(zeta{1}(:,1,idx2)'*pit(1) + zeta{2}(:,1,idx2)'*pit(2) + zeta{3}(:,1,idx2)'*pit(3));
%     end
% end

